
# Load packages
  library(tidyverse)
  library(haven)
  library(broom)
  library(lmtest)
  library(sandwich)
  library(stargazer)
  
# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]
  
  ##########
  ## 2017 ##
  ##########
  
  #### Langfrist-Online-Tracking, Kumulation 2009-2017 ----
  data <- read_dta("ObsStudy_00_data_Germany_GLES_longterm_online_tracking_2013_2017.dta")
  
  # Select 2017 data, waves 35-37
  dt17 <- data[data$sample == 35 | data$sample == 36 | data$sample == 37, ]
  
  # Set values smaller 0 to NA as those values simply categorize reasons for missing values
  dt17[dt17 < 0] <- NA
  
  # Evaluation political parties
  colnames(dt17)[seq(119, 126, 1)] <- c("rat_CDU", "rat_CSU", "rat_SPD", "rat_FDP", "rat_LINKE", "rat_GRUENE", "rat_PIRATEN", "rat_AfD")
  
  # Coalition evaluation
  colnames(dt17)[seq(485, 495, 1)] <- c("rat_coal_CDU_CSU", "rat_coal_CDU_CSU_SPD", "rat_coal_SPD_GRUENE", "rat_coal_CDU_CSU_FDP", "rat_coal_SPD_GRUENE_LINKE", "rat_coal_CDU_CSU_GRUENE", "rat_coal_CDU_CSU_FDP_GRUENE", "rat_coal_SPD", "rat_coal_SPD_FDP_GRUENE", "rat_coal_CDU_CSU_AfD","rat_coal_SPD_FDP")
  
  # Coalition likelihood
  colnames(dt17)[seq(980, 988, 1)] <- c("coallik_CDU_CSU_SPD", "coallik_SPD_GRUENE", "coallik_CDU_CSU_FDP", "coallik_SPD_GRUENE_LINKE", "coallik_CDU_CSU_GRUENE", "coallik_CDU_CSU_GRUENE_minority", "coallik_CDU_CSU_FDP_GRUENE", "coallik_SPD_FDP_GRUENE", "coallik_CDU_CSU_minority")
  
  # Dependent variable
  vote <- data.frame(dt17$k0080bb_1, dt17$k0080bb_2, dt17$k0086bb, dt17$k0090bb_1, dt17$k0090bb_2)

  # Generate party choice variable
  dt17$vote <- dt17$k0080bb_1
  dt17$vote <- ifelse(is.na(dt17$vote), dt17$k0080bb_2, dt17$vote)
  dt17$vote <- ifelse(is.na(dt17$vote), dt17$k0086bb, dt17$vote)
  dt17$vote <- ifelse(is.na(dt17$vote), dt17$k0090bb_1, dt17$vote)
  dt17$vote <- ifelse(is.na(dt17$vote), dt17$k0090bb_2, dt17$vote)

  # Create a "rat_coal_AfD" variable that is equal to rat_AfD, because no other party intends to form a coalition with AfD
  dt17$rat_coal_AfD <- dt17$rat_AfD

  # Assign interpretable labels to vote variable
  dt17$vote[dt17$vote==1] <- "UNION"
  dt17$vote[dt17$vote==4] <- "SPD"
  dt17$vote[dt17$vote==5] <- "FDP"
  dt17$vote[dt17$vote==6] <- "GRUENE"
  dt17$vote[dt17$vote==7] <- "LINKE"
  dt17$vote[dt17$vote==322] <- "AfD"
  
  # Generate ptv columns based on vote variable
  dt17$ptv_UNION <- ifelse(dt17$vote == "UNION", 1, 0)
  dt17$ptv_SPD <- ifelse(dt17$vote == "SPD", 1, 0)
  dt17$ptv_FDP <- ifelse(dt17$vote == "FDP", 1, 0)
  dt17$ptv_GRUENE <- ifelse(dt17$vote == "GRUENE", 1, 0)
  dt17$ptv_LINKE <- ifelse(dt17$vote == "LINKE", 1, 0)
  dt17$ptv_AfD <- ifelse(dt17$vote == "AfD", 1, 0)
  
  Z <- dt17 %>% select("rat_coal_CDU_CSU_SPD","rat_coal_CDU_CSU_FDP", "rat_coal_CDU_CSU_FDP_GRUENE", "rat_coal_CDU_CSU_GRUENE", "rat_coal_SPD_GRUENE_LINKE", "rat_coal_AfD") %>% 
    mutate_all(as.numeric) %>% 
    mutate_all(funs(scales::rescale(.,to = c(0, 1), from = c(1,11)))) 
  
  ##### CDU/CSU #####
  # Coalition evaluation
  rat_coal_CDU_CSU <- Z %>% select(-c(rat_coal_SPD_GRUENE_LINKE, rat_coal_AfD)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood
  coal_lik_CDU_CSU <- dt17 %>% select(coallik_CDU_CSU_SPD, coallik_CDU_CSU_FDP, coallik_CDU_CSU_FDP_GRUENE, coallik_CDU_CSU_GRUENE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coallik_CDU_CSU_SPD) + exp(coallik_CDU_CSU_FDP) + exp(coallik_CDU_CSU_FDP_GRUENE) + exp(coallik_CDU_CSU_GRUENE),
           coallik_CDU_CSU_SPD = exp(coallik_CDU_CSU_SPD)/sum_exp,
           coallik_CDU_CSU_FDP = exp(coallik_CDU_CSU_FDP)/sum_exp,
           coallik_CDU_CSU_FDP_GRUENE = exp(coallik_CDU_CSU_FDP_GRUENE)/sum_exp,
           coallik_CDU_CSU_GRUENE = exp(coallik_CDU_CSU_GRUENE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### SPD #####
  # Coalition evaluation
  rat_coal_SPD <- Z %>% select(-c(rat_coal_CDU_CSU_FDP, rat_coal_CDU_CSU_FDP_GRUENE, rat_coal_CDU_CSU_GRUENE, rat_coal_AfD)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood
  coal_lik_SPD <- dt17 %>% select(coallik_CDU_CSU_SPD, coallik_SPD_GRUENE_LINKE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coallik_CDU_CSU_SPD) + exp(coallik_SPD_GRUENE_LINKE),
           coallik_CDU_CSU_SPD = exp(coallik_CDU_CSU_SPD)/sum_exp,
           coallik_SPD_GRUENE_LINKE = exp(coallik_SPD_GRUENE_LINKE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### FDP #####
  # Coalition evaluation
  rat_coal_FDP <- Z %>% select(-c(rat_coal_CDU_CSU_SPD, rat_coal_CDU_CSU_GRUENE, rat_coal_SPD_GRUENE_LINKE, rat_coal_AfD)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  
  # Coalition likelihood
  coal_lik_FDP <- dt17 %>% select(coallik_CDU_CSU_FDP, coallik_CDU_CSU_FDP_GRUENE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coallik_CDU_CSU_FDP) + exp(coallik_CDU_CSU_FDP_GRUENE),
           coallik_CDU_CSU_FDP = exp(coallik_CDU_CSU_FDP)/sum_exp,
           coallik_CDU_CSU_FDP_GRUENE = exp(coallik_CDU_CSU_FDP_GRUENE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### GRUENE #####
  # Coalition evaluation
  rat_coal_GRUENE <- Z %>% select(-c(rat_coal_CDU_CSU_SPD, rat_coal_CDU_CSU_FDP, rat_coal_AfD)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood
  coal_lik_GRUENE <- dt17 %>% select(coallik_CDU_CSU_FDP_GRUENE, coallik_CDU_CSU_GRUENE, coallik_SPD_GRUENE_LINKE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp =  exp(coallik_CDU_CSU_FDP_GRUENE) + exp(coallik_CDU_CSU_GRUENE) + exp(coallik_SPD_GRUENE_LINKE),
           coallik_CDU_CSU_FDP_GRUENE = exp(coallik_CDU_CSU_FDP_GRUENE)/sum_exp,
           coallik_CDU_CSU_GRUENE = exp(coallik_CDU_CSU_GRUENE)/sum_exp,
           coallik_SPD_GRUENE_LINKE = exp(coallik_SPD_GRUENE_LINKE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### LINKE #####
  # Coalition evaluation
  rat_coal_LINKE <- Z %>% select(rat_coal_SPD_GRUENE_LINKE) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood, has to be 1 as only one coalition option available
  coal_lik_LINKE <- dt17 %>% select(coallik_SPD_GRUENE_LINKE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp =  exp(coallik_SPD_GRUENE_LINKE),
           coallik_SPD_GRUENE_LINKE = exp(coallik_SPD_GRUENE_LINKE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  ##### AfD #####
  # Coalition evaluation
  rat_coal_AfD <- Z %>% select(rat_coal_AfD) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood is 0 as no other party would enter coalition with AfD
  coal_lik_AfD <- matrix(rep(0, times = nrow(rat_coal_AfD)))
  
  
  # Mean-Variance Model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(coal_lik_CDU_CSU)
  M_CDU_CSU <- sapply(1:N, function(i) EV(coal_lik_CDU_CSU[i,], rat_coal_CDU_CSU[i,]))
  V_CDU_CSU <- sapply(1:N, function(i) Vari(coal_lik_CDU_CSU[i,], rat_coal_CDU_CSU[i,]))
  
  N <- nrow(coal_lik_SPD)
  M_SPD <- sapply(1:N, function(i) EV(coal_lik_SPD[i,], rat_coal_SPD[i,]))
  V_SPD <- sapply(1:N, function(i) Vari(coal_lik_SPD[i,], rat_coal_SPD[i,]))
  
  N <- nrow(coal_lik_FDP)
  M_FDP <- sapply(1:N, function(i) EV(coal_lik_FDP[i,], rat_coal_FDP[i,]))
  V_FDP <- sapply(1:N, function(i) Vari(coal_lik_FDP[i,], rat_coal_FDP[i,]))
  
  N <- nrow(coal_lik_GRUENE)
  M_GRUENE <- sapply(1:N, function(i) EV(coal_lik_GRUENE[i,], rat_coal_GRUENE[i,]))
  V_GRUENE <- sapply(1:N, function(i) Vari(coal_lik_GRUENE[i,], rat_coal_GRUENE[i,]))
  
  N <- nrow(coal_lik_LINKE)
  M_LINKE <- sapply(1:N, function(i) EV(coal_lik_LINKE[i,], rat_coal_LINKE[i,]))
  V_LINKE <- sapply(1:N, function(i) Vari(coal_lik_LINKE[i,], rat_coal_LINKE[i,])) # is 0
  
  N <- nrow(coal_lik_AfD)
  M_AfD <- sapply(1:N, function(i) EV(coal_lik_AfD[i,], rat_coal_AfD[i,])) # mean 0
  V_AfD <- sapply(1:N, function(i) Vari(coal_lik_AfD[i,], rat_coal_AfD[i,])) # variance 0
  
  # Generate a gender dummy "male"
  dt17$male <- ifelse(dt17$k0010 == 2, 0, dt17$k0010)

  # Age
  dt17$age <- 2017 - dt17$k0020
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt17$pid_UNION <- ifelse(dt17$k0490b==1|dt17$k0490b==2|dt17$k0490b==3, 1, 0)
  dt17$pid_SPD <- ifelse(dt17$k0490b==4, 1, 0)
  dt17$pid_FDP <- ifelse(dt17$k0490b==5, 1, 0)
  dt17$pid_GRUENE <- ifelse(dt17$k0490b==6, 1, 0)
  dt17$pid_LINKE <- ifelse(dt17$k0490b==7, 1, 0)
  dt17$pid_AfD <- ifelse(dt17$k0490b==322, 1, 0)
  
  d_gles_ot_2017 <- data.frame(
    "ptv" = dt17$vote,
    "ptv_UNION" = dt17$ptv_UNION,
    "ptv_SPD" = dt17$ptv_SPD,
    "ptv_FDP" = dt17$ptv_FDP,
    "ptv_GRUENE" = dt17$ptv_GRUENE,
    "ptv_LINKE" = dt17$ptv_LINKE,
    "ptv_AfD" = dt17$ptv_AfD,
    "rat.UNION" = (dt17$rat_CDU + dt17$rat_CSU)/2,
    "rat.SPD" = dt17$rat_SPD,
    "rat.FDP" = dt17$rat_FDP,
    "rat.GRUENE" = dt17$rat_GRUENE,
    "rat.LINKE" = dt17$rat_LINKE,
    "rat.AfD" = dt17$rat_AfD,
    "pid_UNION" = dt17$pid_UNION,
    "pid_SPD" = dt17$pid_SPD,
    "pid_FDP" = dt17$pid_FDP,
    "pid_GRUENE" = dt17$pid_GRUENE,
    "pid_LINKE" = dt17$pid_LINKE,
    "pid_AfD" = dt17$pid_AfD,
    "lotterymean.UNION" = M_CDU_CSU,
    "lotteryvariance.UNION" = V_CDU_CSU,
    "lotterymean.SPD" = M_SPD,
    "lotteryvariance.SPD" = V_SPD,
    "lotterymean.FDP" = M_FDP,
    "lotteryvariance.FDP" = V_FDP,
    "lotterymean.GRUENE" = M_GRUENE,
    "lotteryvariance.GRUENE" = V_GRUENE,
    "lotterymean.LINKE" = M_LINKE,
    "lotteryvariance.LINKE" = V_LINKE,
    "lotterymean.AfD" = M_AfD,
    "lotteryvariance.AfD" = V_AfD,
    "sex" = dt17$male,
    "age" = dt17$age,
    "edu" = dt17$k0030_2)
  
  d_gles_ot_2017$lotteryvariance.UNION <- scales::rescale(d_gles_ot_2017$lotteryvariance.UNION, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2017$lotteryvariance.SPD <- scales::rescale(d_gles_ot_2017$lotteryvariance.SPD, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2017$lotteryvariance.FDP <- scales::rescale(d_gles_ot_2017$lotteryvariance.FDP, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2017$lotteryvariance.GRUENE <- scales::rescale(d_gles_ot_2017$lotteryvariance.GRUENE, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2017$lotteryvariance.LINKE <- scales::rescale(d_gles_ot_2017$lotteryvariance.LINKE, to = c(0, 1), from = c(0,0.25))
  
# Saving model results
  summary(m1.gles.ot.2017 <- lm(as.formula(paste("ptv_UNION ~ lotteryvariance.UNION + lotterymean.UNION + rat.UNION  + sex + as.factor(edu) + age", 
                                                 if (robustnesscheck_pid) "+ pid_UNION" else "")), d_gles_ot_2017))
  summary(m2.gles.ot.2017 <- lm(as.formula(paste("ptv_SPD ~ lotteryvariance.SPD + lotterymean.SPD + rat.SPD + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_SPD" else "")), d_gles_ot_2017))
  summary(m3.gles.ot.2017 <- lm(as.formula(paste("ptv_FDP ~ lotteryvariance.FDP + lotterymean.FDP + rat.FDP + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_FDP" else "")), d_gles_ot_2017))
  summary(m4.gles.ot.2017 <- lm(as.formula(paste("ptv_GRUENE ~ lotteryvariance.GRUENE + lotterymean.GRUENE + rat.GRUENE + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_GRUENE" else "")), d_gles_ot_2017))
  summary(m5.gles.ot.2017 <- lm(as.formula(paste("ptv_LINKE ~ lotteryvariance.LINKE + lotterymean.LINKE + rat.LINKE + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_LINKE" else "")), d_gles_ot_2017))
  
  GER_2017_GLESLOT_UNION_Variance_Estimate <- tidy(m1.gles.ot.2017) %>% filter(term=="lotteryvariance.UNION") %>% select("estimate")
  GER_2017_GLESLOT_UNION_Variance_SE <- se_robust(m1.gles.ot.2017)["lotteryvariance.UNION"] #tidy(m1.gles.ot.2017) %>% filter(term=="lotteryvariance.UNION") %>% select("std.error")
  GER_2017_GLESLOT_UNION_Mean_Estimate <- tidy(m1.gles.ot.2017) %>% filter(term=="lotterymean.UNION") %>% select("estimate")
  GER_2017_GLESLOT_UNION_Mean_SE <- se_robust(m1.gles.ot.2017)["lotterymean.UNION"] #tidy(m1.gles.ot.2017) %>% filter(term=="lotterymean.UNION") %>% select("std.error")
  
  GER_2017_GLESLOT_SPD_Variance_Estimate <- tidy(m2.gles.ot.2017) %>% filter(term=="lotteryvariance.SPD") %>% select("estimate")
  GER_2017_GLESLOT_SPD_Variance_SE <- se_robust(m2.gles.ot.2017)["lotteryvariance.SPD"] #tidy(m2.gles.ot.2017) %>% filter(term=="lotteryvariance.SPD") %>% select("std.error")
  GER_2017_GLESLOT_SPD_Mean_Estimate <- tidy(m2.gles.ot.2017) %>% filter(term=="lotterymean.SPD") %>% select("estimate")
  GER_2017_GLESLOT_SPD_Mean_SE <- se_robust(m2.gles.ot.2017)["lotterymean.SPD"] #tidy(m2.gles.ot.2017) %>% filter(term=="lotterymean.SPD") %>% select("std.error")
  
  GER_2017_GLESLOT_FDP_Variance_Estimate <- tidy(m3.gles.ot.2017) %>% filter(term=="lotteryvariance.FDP") %>% select("estimate")
  GER_2017_GLESLOT_FDP_Variance_SE <- se_robust(m3.gles.ot.2017)["lotteryvariance.FDP"] #tidy(m3.gles.ot.2017) %>% filter(term=="lotteryvariance.FDP") %>% select("std.error")
  GER_2017_GLESLOT_FDP_Mean_Estimate <- tidy(m3.gles.ot.2017) %>% filter(term=="lotterymean.FDP") %>% select("estimate")
  GER_2017_GLESLOT_FDP_Mean_SE <- se_robust(m3.gles.ot.2017)["lotterymean.FDP"] #tidy(m3.gles.ot.2017) %>% filter(term=="lotterymean.FDP") %>% select("std.error")
  
  GER_2017_GLESLOT_GR_Variance_Estimate <- tidy(m4.gles.ot.2017) %>% filter(term=="lotteryvariance.GRUENE") %>% select("estimate")
  GER_2017_GLESLOT_GR_Variance_SE <- se_robust(m4.gles.ot.2017)["lotteryvariance.GRUENE"] #tidy(m4.gles.ot.2017) %>% filter(term=="lotteryvariance.GRUENE") %>% select("std.error")
  GER_2017_GLESLOT_GR_Mean_Estimate <- tidy(m4.gles.ot.2017) %>% filter(term=="lotterymean.GRUENE") %>% select("estimate")
  GER_2017_GLESLOT_GR_Mean_SE <- se_robust(m4.gles.ot.2017)["lotterymean.GRUENE"] #tidy(m4.gles.ot.2017) %>% filter(term=="lotterymean.GRUENE") %>% select("std.error")
  
  # Harmonize names of the IVs of interest in the models
  names(m1.gles.ot.2017$coefficients)[names(m1.gles.ot.2017$coefficients) == "lotteryvariance.UNION"] <- "Government Lottery Variance"
  names(m1.gles.ot.2017$coefficients)[names(m1.gles.ot.2017$coefficients) == "lotterymean.UNION"] <- "Government Lottery Mean"
  names(m2.gles.ot.2017$coefficients)[names(m2.gles.ot.2017$coefficients) == "lotteryvariance.SPD"] <- "Government Lottery Variance"
  names(m2.gles.ot.2017$coefficients)[names(m2.gles.ot.2017$coefficients) == "lotterymean.SPD"] <- "Government Lottery Mean"
  names(m3.gles.ot.2017$coefficients)[names(m3.gles.ot.2017$coefficients) == "lotteryvariance.FDP"] <- "Government Lottery Variance"
  names(m3.gles.ot.2017$coefficients)[names(m3.gles.ot.2017$coefficients) == "lotterymean.FDP"] <- "Government Lottery Mean"
  names(m4.gles.ot.2017$coefficients)[names(m4.gles.ot.2017$coefficients) == "lotteryvariance.GRUENE"] <- "Government Lottery Variance"
  names(m4.gles.ot.2017$coefficients)[names(m4.gles.ot.2017$coefficients) == "lotterymean.GRUENE"] <- "Government Lottery Mean"
  names(m5.gles.ot.2017$coefficients)[names(m5.gles.ot.2017$coefficients) == "lotteryvariance.LINKE"] <- "Government Lottery Variance"
  names(m5.gles.ot.2017$coefficients)[names(m5.gles.ot.2017$coefficients) == "lotterymean.LINKE"] <- "Government Lottery Mean"
  
  if(!robustnesscheck_pid){
    stargazer(m1.gles.ot.2017, m2.gles.ot.2017, m3.gles.ot.2017, m4.gles.ot.2017, title="Linear regressions of vote choice on perceived government lottery variance and mean (GLES Longterm-Online-Tracking 2017).", 
              dep.var.labels=c("Union", "SPD", "FDP", "Greens", "Left"), no.space=T,
              omit.stat=c("LL","ser","f"),
              se = lapply(list(m1.gles.ot.2017, m2.gles.ot.2017, m3.gles.ot.2017, m4.gles.ot.2017), se_robust),
              p = lapply(list(m1.gles.ot.2017, m2.gles.ot.2017, m3.gles.ot.2017, m4.gles.ot.2017), p_robust),
              type = "text",
              out = "TableSM9.tex")
    }

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  ##########
  ## 2013 ##
  ##########
  
  #### Langfrist-Online-Tracking, Kumulation 2009-2017
  data <- read_dta("ObsStudy_00_data_Germany_GLES_longterm_online_tracking_2013_2017.dta")
  
  # Select 2013 data, waves 19, 20 and 21
  dt13 <- data[data$sample == 19 | data$sample == 20 |data$sample == 21, ]
  
  # Set values smaller 0 to NA as those values simply categorize reasons for missing values
  dt13[dt13 < 0] <- NA
  
  # Evaluation political parties
  colnames(dt13)[seq(119, 126, 1)] <- c("rat_CDU", "rat_CSU", "rat_SPD", "rat_FDP", "rat_LINKE", "rat_GRUENE", "rat_PIRATEN", "rat_AfD")
  
  # Coalition evaluation
  colnames(dt13)[seq(485, 493, 1)] <- c("rat_coal_CDU_CSU", "rat_coal_CDU_CSU_SPD", "rat_coal_SPD_GRUENE", "rat_coal_CDU_CSU_FDP", "rat_coal_SPD_GRUENE_LINKE", "rat_coal_CDU_CSU_GRUENE", "rat_coal_CDU_CSU_FDP_GRUENE", "rat_coal_SPD", "rat_coal_SPD_FDP_GRUENE")

  # Coalition likelihood
  colnames(dt13)[seq(980, 988, 1)] <- c("coallik_CDU_CSU_SPD", "coallik_SPD_GRUENE", "coallik_CDU_CSU_FDP", "coallik_SPD_GRUENE_LINKE", "coallik_CDU_CSU_GRUENE", "coallik_CDU_CSU_GRUENE_minority", "coallik_CDU_CSU_FDP_GRUENE", "coallik_SPD_FDP_GRUENE", "coallik_CDU_CSU_minority")
  
  summary(dt13$coallik_CDU_CSU_SPD)
  summary(dt13$coallik_SPD_GRUENE)
  summary(dt13$coallik_CDU_CSU_FDP)
  summary(dt13$coallik_SPD_GRUENE_LINKE)
  summary(dt13$coallik_CDU_CSU_GRUENE)
  summary(dt13$coallik_CDU_CSU_GRUENE_minority) # not asked
  summary(dt13$coallik_CDU_CSU_FDP_GRUENE)
  summary(dt13$coallik_SPD_FDP_GRUENE)
  summary(dt13$coallik_CDU_CSU_minority) # not asked
  
  # Dependent variable
  vote <- data.frame(dt13$k0080bb_1, dt13$k0080bb_2, dt13$k0086bb, dt13$k0090bb_1, dt13$k0090bb_2)

  dt13$vote <- dt13$k0080bb_1
  dt13$vote <- ifelse(is.na(dt13$vote), dt13$k0080bb_2, dt13$vote)
  dt13$vote <- ifelse(is.na(dt13$vote), dt13$k0086bb, dt13$vote)
  dt13$vote <- ifelse(is.na(dt13$vote), dt13$k0090bb_1, dt13$vote)
  dt13$vote <- ifelse(is.na(dt13$vote), dt13$k0090bb_2, dt13$vote)

  # Assign interpretable labels to vote variable
  dt13$vote[dt13$vote==1] <- "UNION"
  dt13$vote[dt13$vote==4] <- "SPD"
  dt13$vote[dt13$vote==5] <- "FDP"
  dt13$vote[dt13$vote==6] <- "GRUENE"
  dt13$vote[dt13$vote==7] <- "LINKE"
  
  ## Generate ptv columns based on vote variable
  dt13$ptv_UNION <- ifelse(dt13$vote == "UNION", 1, 0)
  dt13$ptv_SPD <- ifelse(dt13$vote == "SPD", 1, 0)
  dt13$ptv_FDP <- ifelse(dt13$vote == "FDP", 1, 0)
  dt13$ptv_GRUENE <- ifelse(dt13$vote == "GRUENE", 1, 0)
  dt13$ptv_LINKE <- ifelse(dt13$vote == "LINKE", 1, 0)
  
  Z <- dt13 %>% select("rat_coal_CDU_CSU_SPD", "rat_coal_SPD_GRUENE" ,"rat_coal_CDU_CSU_FDP", "rat_coal_SPD_GRUENE_LINKE", "rat_coal_CDU_CSU_GRUENE", "rat_coal_CDU_CSU_FDP_GRUENE", "rat_coal_SPD_FDP_GRUENE") %>% 
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1), from = c(1,11)))) 
  
  
  ##### CDU/CSU #####
  # Coalition evaluation
  rat_coal_CDU_CSU <- Z %>% select(-c(rat_coal_SPD_GRUENE, rat_coal_SPD_GRUENE_LINKE, rat_coal_SPD_FDP_GRUENE)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood
  coal_lik_CDU_CSU <- dt13 %>% select(coallik_CDU_CSU_SPD, coallik_CDU_CSU_FDP, coallik_CDU_CSU_FDP_GRUENE, coallik_CDU_CSU_GRUENE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coallik_CDU_CSU_SPD) + exp(coallik_CDU_CSU_FDP) + exp(coallik_CDU_CSU_FDP_GRUENE) + exp(coallik_CDU_CSU_GRUENE),
           coallik_CDU_CSU_SPD = exp(coallik_CDU_CSU_SPD)/sum_exp,
           coallik_CDU_CSU_FDP = exp(coallik_CDU_CSU_FDP)/sum_exp,
           coallik_CDU_CSU_FDP_GRUENE = exp(coallik_CDU_CSU_FDP_GRUENE)/sum_exp,
           coallik_CDU_CSU_GRUENE = exp(coallik_CDU_CSU_GRUENE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### SPD #####
  # Coalition evaluation
  rat_coal_SPD <- Z %>% select(-c(rat_coal_CDU_CSU_FDP, rat_coal_CDU_CSU_FDP_GRUENE, rat_coal_CDU_CSU_GRUENE)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood
  coal_lik_SPD <- dt13 %>% select(coallik_CDU_CSU_SPD, coallik_SPD_GRUENE, coallik_SPD_FDP_GRUENE, coallik_SPD_GRUENE_LINKE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coallik_CDU_CSU_SPD) + exp(coallik_SPD_GRUENE) + exp(coallik_SPD_FDP_GRUENE) + exp(coallik_SPD_GRUENE_LINKE),
           coallik_CDU_CSU_SPD = exp(coallik_CDU_CSU_SPD)/sum_exp,
           coallik_SPD_GRUENE = exp(coallik_SPD_GRUENE)/sum_exp,
           coallik_SPD_FDP_GRUENE = exp(coallik_SPD_FDP_GRUENE)/sum_exp,
           coallik_SPD_GRUENE_LINKE = exp(coallik_SPD_GRUENE_LINKE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### FDP #####
  # Coalition evaluation
  rat_coal_FDP <- Z %>% select(-c(rat_coal_CDU_CSU_SPD, rat_coal_SPD_GRUENE, rat_coal_CDU_CSU_GRUENE, rat_coal_SPD_GRUENE_LINKE)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  
  # Coalition likelihood
  coal_lik_FDP <- dt13 %>% select(coallik_CDU_CSU_FDP, coallik_SPD_FDP_GRUENE, coallik_CDU_CSU_FDP_GRUENE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coallik_CDU_CSU_FDP) + exp(coallik_SPD_FDP_GRUENE) + exp(coallik_CDU_CSU_FDP_GRUENE),
           coallik_CDU_CSU_FDP = exp(coallik_CDU_CSU_FDP)/sum_exp,
           coallik_SPD_FDP_GRUENE = exp(coallik_SPD_FDP_GRUENE)/sum_exp,
           coallik_CDU_CSU_FDP_GRUENE = exp(coallik_CDU_CSU_FDP_GRUENE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### GRUENE #####
  # Coalition evaluation
  rat_coal_GRUENE <- Z %>% select(-c(rat_coal_CDU_CSU_SPD, rat_coal_CDU_CSU_FDP)) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood
  coal_lik_GRUENE <- dt13 %>% select(coallik_SPD_GRUENE, coallik_CDU_CSU_FDP_GRUENE, coallik_CDU_CSU_GRUENE, coallik_SPD_FDP_GRUENE, coallik_SPD_GRUENE_LINKE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp =  exp(coallik_SPD_GRUENE) + exp(coallik_CDU_CSU_FDP_GRUENE) + exp(coallik_CDU_CSU_GRUENE) + exp(coallik_SPD_FDP_GRUENE) + exp(coallik_SPD_GRUENE_LINKE),
           coallik_SPD_GRUENE = exp(coallik_SPD_GRUENE)/sum_exp,
           coallik_CDU_CSU_FDP_GRUENE = exp(coallik_CDU_CSU_FDP_GRUENE)/sum_exp,
           coallik_CDU_CSU_GRUENE = exp(coallik_CDU_CSU_GRUENE)/sum_exp,
           coallik_SPD_FDP_GRUENE = exp(coallik_SPD_FDP_GRUENE)/sum_exp,
           coallik_SPD_GRUENE_LINKE = exp(coallik_SPD_GRUENE_LINKE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  ##### LINKE #####
  # Coalition evaluation
  rat_coal_LINKE <- Z %>% select(rat_coal_SPD_GRUENE_LINKE) %>%
    mutate_all(as.numeric) %>%
    as.matrix()
  
  # Coalition likelihood, has to be 1 as only one coalition option available
  coal_lik_LINKE <- dt13 %>% select(coallik_SPD_GRUENE_LINKE)  %>%
    mutate_all(as.numeric) %>%
    mutate(sum_exp =  exp(coallik_SPD_GRUENE_LINKE),
           coallik_SPD_GRUENE_LINKE = exp(coallik_SPD_GRUENE_LINKE)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  # Mean-Variance model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  
  N <- nrow(coal_lik_CDU_CSU)
  M_CDU_CSU <- sapply(1:N, function(i) EV(coal_lik_CDU_CSU[i,], rat_coal_CDU_CSU[i,]))
  V_CDU_CSU <- sapply(1:N, function(i) Vari(coal_lik_CDU_CSU[i,], rat_coal_CDU_CSU[i,]))
  
  N <- nrow(coal_lik_SPD)
  M_SPD <- sapply(1:N, function(i) EV(coal_lik_SPD[i,], rat_coal_SPD[i,]))
  V_SPD <- sapply(1:N, function(i) Vari(coal_lik_SPD[i,], rat_coal_SPD[i,]))
  
  N <- nrow(coal_lik_FDP)
  M_FDP <- sapply(1:N, function(i) EV(coal_lik_FDP[i,], rat_coal_FDP[i,]))
  V_FDP <- sapply(1:N, function(i) Vari(coal_lik_FDP[i,], rat_coal_FDP[i,]))
  
  N <- nrow(coal_lik_GRUENE)
  M_GRUENE <- sapply(1:N, function(i) EV(coal_lik_GRUENE[i,], rat_coal_GRUENE[i,]))
  V_GRUENE <- sapply(1:N, function(i) Vari(coal_lik_GRUENE[i,], rat_coal_GRUENE[i,]))
  
  N <- nrow(coal_lik_LINKE)
  M_LINKE <- sapply(1:N, function(i) EV(coal_lik_LINKE[i,], rat_coal_LINKE[i,]))
  V_LINKE <- sapply(1:N, function(i) Vari(coal_lik_LINKE[i,], rat_coal_LINKE[i,]))
  
  # Generate a gender dummy "male"
  dt13$male <- ifelse(dt13$k0010 == 2, 0, dt13$k0010)

  # Age
  dt13$age <- 2013 - dt13$k0020
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt13$pid_UNION <- ifelse(dt13$k0490b==1|dt13$k0490b==2|dt13$k0490b==3, 1, 0)
  dt13$pid_SPD <- ifelse(dt13$k0490b==4, 1, 0)
  dt13$pid_FDP <- ifelse(dt13$k0490b==5, 1, 0)
  dt13$pid_GRUENE <- ifelse(dt13$k0490b==6, 1, 0)
  dt13$pid_LINKE <- ifelse(dt13$k0490b==7, 1, 0)
  dt13$pid_AfD <- ifelse(dt13$k0490b==322, 1, 0)

  d_gles_ot_2013 <- data.frame(
    "ptv" = dt13$vote,
    "ptv_UNION" = dt13$ptv_UNION,
    "ptv_SPD" = dt13$ptv_SPD,
    "ptv_FDP" = dt13$ptv_FDP,
    "ptv_GRUENE" = dt13$ptv_GRUENE,
    "ptv_LINKE" = dt13$ptv_LINKE,
    "rat.UNION" = (dt13$rat_CDU + dt13$rat_CSU)/2,
    "rat.SPD" = dt13$rat_SPD,
    "rat.FDP" = dt13$rat_FDP,
    "rat.GRUENE" = dt13$rat_GRUENE,
    "rat.LINKE" = dt13$rat_LINKE,
    "pid_UNION" = dt13$pid_UNION,
    "pid_SPD" = dt13$pid_SPD,
    "pid_FDP" = dt13$pid_FDP,
    "pid_GRUENE" = dt13$pid_GRUENE,
    "pid_LINKE" = dt13$pid_LINKE,
    "pid_AfD" = dt13$pid_AfD,
    "lotterymean.UNION" = M_CDU_CSU,
    "lotteryvariance.UNION" = V_CDU_CSU,
    "lotterymean.SPD" = M_SPD,
    "lotteryvariance.SPD" = V_SPD,
    "lotterymean.FDP" = M_FDP,
    "lotteryvariance.FDP" = V_FDP,
    "lotterymean.GRUENE" = M_GRUENE,
    "lotteryvariance.GRUENE" = V_GRUENE,
    "lotterymean.LINKE" = M_LINKE,
    "lotteryvariance.LINKE" = V_LINKE,
    "sex" = dt13$male,
    "age" = dt13$age,
    "edu" = dt13$k0030_2
    )
  
  d_gles_ot_2013$lotteryvariance.UNION <- scales::rescale(d_gles_ot_2013$lotteryvariance.UNION, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2013$lotteryvariance.SPD <- scales::rescale(d_gles_ot_2013$lotteryvariance.SPD, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2013$lotteryvariance.FDP <- scales::rescale(d_gles_ot_2013$lotteryvariance.FDP, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2013$lotteryvariance.GRUENE <- scales::rescale(d_gles_ot_2013$lotteryvariance.GRUENE, to = c(0, 1), from = c(0,0.25))
  d_gles_ot_2013$lotteryvariance.LINKE <- scales::rescale(d_gles_ot_2013$lotteryvariance.LINKE, to = c(0, 1), from = c(0,0.25))
  
# Saving model results
  summary(m1.gles.ot.2013 <- lm(as.formula(paste("ptv_UNION ~ lotteryvariance.UNION + lotterymean.UNION + rat.UNION  + sex + as.factor(edu) + age", 
                                                 if (robustnesscheck_pid) "+ pid_UNION" else "")), d_gles_ot_2013))
  summary(m2.gles.ot.2013 <- lm(as.formula(paste("ptv_SPD ~ lotteryvariance.SPD + lotterymean.SPD + rat.SPD + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_SPD" else "")), d_gles_ot_2013))
  summary(m3.gles.ot.2013 <- lm(as.formula(paste("ptv_FDP ~ lotteryvariance.FDP + lotterymean.FDP + rat.FDP + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_FDP" else "")), d_gles_ot_2013))
  summary(m4.gles.ot.2013 <- lm(as.formula(paste("ptv_GRUENE ~ lotteryvariance.GRUENE + lotterymean.GRUENE + rat.GRUENE + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_GRUENE" else "")), d_gles_ot_2013))
  summary(m5.gles.ot.2013 <- lm(as.formula(paste("ptv_LINKE ~ lotteryvariance.LINKE + lotterymean.LINKE + rat.LINKE + sex + as.factor(edu) + age",
                                                 if (robustnesscheck_pid) "+ pid_LINKE" else "")), d_gles_ot_2013))
  
  GER_2013_GLESLOT_UNION_Variance_Estimate <- tidy(m1.gles.ot.2013) %>% filter(term=="lotteryvariance.UNION") %>% select("estimate")
  GER_2013_GLESLOT_UNION_Variance_SE <- se_robust(m1.gles.ot.2013)["lotteryvariance.UNION"] #tidy(m1.gles.ot.2013) %>% filter(term=="lotteryvariance.UNION") %>% select("std.error")
  GER_2013_GLESLOT_UNION_Mean_Estimate <- tidy(m1.gles.ot.2013) %>% filter(term=="lotterymean.UNION") %>% select("estimate")
  GER_2013_GLESLOT_UNION_Mean_SE <- se_robust(m1.gles.ot.2013)["lotterymean.UNION"] #tidy(m1.gles.ot.2013) %>% filter(term=="lotterymean.UNION") %>% select("std.error")
  
  GER_2013_GLESLOT_SPD_Variance_Estimate <- tidy(m2.gles.ot.2013) %>% filter(term=="lotteryvariance.SPD") %>% select("estimate")
  GER_2013_GLESLOT_SPD_Variance_SE <- se_robust(m2.gles.ot.2013)["lotteryvariance.SPD"] #tidy(m2.gles.ot.2013) %>% filter(term=="lotteryvariance.SPD") %>% select("std.error")
  GER_2013_GLESLOT_SPD_Mean_Estimate <- tidy(m2.gles.ot.2013) %>% filter(term=="lotterymean.SPD") %>% select("estimate")
  GER_2013_GLESLOT_SPD_Mean_SE <- se_robust(m2.gles.ot.2013)["lotterymean.SPD"] #tidy(m2.gles.ot.2013) %>% filter(term=="lotterymean.SPD") %>% select("std.error")
  
  GER_2013_GLESLOT_FDP_Variance_Estimate <- tidy(m3.gles.ot.2013) %>% filter(term=="lotteryvariance.FDP") %>% select("estimate")
  GER_2013_GLESLOT_FDP_Variance_SE <- se_robust(m3.gles.ot.2013)["lotteryvariance.FDP"] #tidy(m3.gles.ot.2013) %>% filter(term=="lotteryvariance.FDP") %>% select("std.error")
  GER_2013_GLESLOT_FDP_Mean_Estimate <- tidy(m3.gles.ot.2013) %>% filter(term=="lotterymean.FDP") %>% select("estimate")
  GER_2013_GLESLOT_FDP_Mean_SE <- se_robust(m3.gles.ot.2013)["lotterymean.FDP"] #tidy(m3.gles.ot.2013) %>% filter(term=="lotterymean.FDP") %>% select("std.error")
  
  GER_2013_GLESLOT_GR_Variance_Estimate <- tidy(m4.gles.ot.2013) %>% filter(term=="lotteryvariance.GRUENE") %>% select("estimate")
  GER_2013_GLESLOT_GR_Variance_SE <- se_robust(m4.gles.ot.2013)["lotteryvariance.GRUENE"] #tidy(m4.gles.ot.2013) %>% filter(term=="lotteryvariance.GRUENE") %>% select("std.error")
  GER_2013_GLESLOT_GR_Mean_Estimate <- tidy(m4.gles.ot.2013) %>% filter(term=="lotterymean.GRUENE") %>% select("estimate")
  GER_2013_GLESLOT_GR_Mean_SE <- se_robust(m4.gles.ot.2013)["lotterymean.GRUENE"] #tidy(m4.gles.ot.2013) %>% filter(term=="lotterymean.GRUENE") %>% select("std.error")
  
  # Harmonize names of the IVs of interest in the models
  names(m1.gles.ot.2013$coefficients)[names(m1.gles.ot.2013$coefficients) == "lotteryvariance.UNION"] <- "Government Lottery Variance"
  names(m1.gles.ot.2013$coefficients)[names(m1.gles.ot.2013$coefficients) == "lotterymean.UNION"] <- "Government Lottery Mean"
  names(m2.gles.ot.2013$coefficients)[names(m2.gles.ot.2013$coefficients) == "lotteryvariance.SPD"] <- "Government Lottery Variance"
  names(m2.gles.ot.2013$coefficients)[names(m2.gles.ot.2013$coefficients) == "lotterymean.SPD"] <- "Government Lottery Mean"
  names(m3.gles.ot.2013$coefficients)[names(m3.gles.ot.2013$coefficients) == "lotteryvariance.FDP"] <- "Government Lottery Variance"
  names(m3.gles.ot.2013$coefficients)[names(m3.gles.ot.2013$coefficients) == "lotterymean.FDP"] <- "Government Lottery Mean"
  names(m4.gles.ot.2013$coefficients)[names(m4.gles.ot.2013$coefficients) == "lotteryvariance.GRUENE"] <- "Government Lottery Variance"
  names(m4.gles.ot.2013$coefficients)[names(m4.gles.ot.2013$coefficients) == "lotterymean.GRUENE"] <- "Government Lottery Mean"
  names(m5.gles.ot.2013$coefficients)[names(m5.gles.ot.2013$coefficients) == "lotteryvariance.LINKE"] <- "Government Lottery Variance"
  names(m5.gles.ot.2013$coefficients)[names(m5.gles.ot.2013$coefficients) == "lotterymean.LINKE"] <- "Government Lottery Mean"

  
  if(!robustnesscheck_pid){
  stargazer(m1.gles.ot.2013, m2.gles.ot.2013, m3.gles.ot.2013, m4.gles.ot.2013, title="Linear regressions of vote choice on perceived government lottery variance and mean (GLES Longterm-Online-Tracking 2013).", 
             dep.var.labels=c("Union", "SPD", "FDP", "Greens", "Left"), no.space=T,
            omit.stat=c("LL","ser","f"),
            se = lapply(list(m1.gles.ot.2013, m2.gles.ot.2013, m3.gles.ot.2013, m4.gles.ot.2013), se_robust),
            p = lapply(list(m1.gles.ot.2013, m2.gles.ot.2013, m3.gles.ot.2013, m4.gles.ot.2013), p_robust),
            type = "text",
            out = "TableSM8.tex"
            )
  }
